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Abstract 

Weakly interacting massive particles (WIMPs) are one of the leading candidates for 
Dark Matter. Currently, the most promising method to detect many different WIMP 
candidates is the direct detection of the recoil energy deposited in a low-background lab- 
oratory detector due to elastic WIMP-nucleus scattering. So far the usual procedure has 
been to predict the event rate of direct detection of WIMPs based on some model(s) of 
the galactic halo. The aim of our work is to invert this process. That is, we study what 
future direct detection experiment can teach us about the WIMP halo. As the first step we 
consider a time-averaged recoil spectrum, assuming that no directional information exists. 
We develop a method to construct the (time-averaged) one-dimensional velocity distribu- 
tion function from this spectrum. Moments of this function, such as the mean velocity and 
velocity dispersion of WIMPs, can also be obtained directly from the recoil spectrum. The 
only input needed in addition to a measured recoil spectrum is the mass of the WIMP; no 
information about the scattering cross section or WIMP density is required. 



1 Introduction 



The first indication of the existence of Dark Matter has already been found in the 1930s [1] . By 
now astrophysicists have strong evidence [l]-[4] to believe that a large fraction (more than 80%) 
of the matter in the Universe is dark (i.e., interacts at most very weakly with electromagnetic 
radiation). The dominant component of this cosmological dark matter must be due to some yet 
to be discovered, non-baryonic particles. Weakly interacting massive particles (WIMPs) x are 
one of the leading candidates for Dark Matter. WIMPs are stable particles which arise in several 
extensions of the standard model of electroweak interactions. Typically they are presumed to 
have masses between 10 GeV and a few TeV and interact with ordinary matter only weakly (for 
reviews, see [5]). 

Currently, the most promising method to detect many different WIMP candidates is the 
direct detection of the recoil energy deposited in a low-background laboratory detector by elastic 
scattering of ambient WIMPs on the nuclei in a detector [6]- [8]. The event rate of direct WIMP 
detection depends strongly on the velocity distribution of the incident particles. Usually and 
for simplicity, the local velocity distribution of WIMPs is assumed to be a shifted Maxwell 
distribution, as would arise if the Milky Way halo is isothermal [6]- [9]. However, our halo 
is certainly not a precisely isothermal sphere. Possibilities that have been considered in the 
literature include axisymmetric halo models [10], the so-called secondary infall model of halo 
formation [11], and a possible bulk rotations of the halo of our galaxy [12, 13]. 

If the halo of our galaxy consists of WIMPs, about 10 5 WIMPs should pass through every 
square centimeter of the Earth's (and our!) surface per second (for m x ~ 100 GeV). However, 
the cross section of WIMPs on ordinary materials is very low and makes these interactions quite 
rare [5]. For example, in typical SUSY models with neutralino WIMP, the event rate is about 
10~ 4 ~ 1 event/kg/day and the energy deposited in the detector by a single interaction is about 
1 ~ 100 keV. Typical event rates due to cosmic rays and ambient radioactivity are much larger. 
The annual modulation of the event rate due to the orbital motion of the Earth around the Sun 
has been suggested as a way to discriminate signal from background [7, 8, 14, 13]. Actually, the 
DAMA collaboration has claimed that they have observed this annual modulation of the event 
rate [15]. Note, however, that the annual modulation of the signal is expected to amount only 
to a few percent; this method can therefore only be used once more than one hundred signal 
events have been accumulated. In the meantime, a more promising approach is to reduce the 
background by vetoing events that do not look like nuclear recoil. This method is, e.g., being 
used by the CDMS [16], CRESST [17] and EDELWEISS [18] collaborations. The presently best 
null result, from CDMS [19], contradicts the DAMA claim for standard halo models. 1 

So far most theoretical analyses of direct WIMP detection have predicted the detection rate 
for a given (class of) WIMP(s), based on a specific model of the galactic halo. The goal of this 
paper is to invert this process. That is, we wish to study, as model-independently as possible, 
what future direct detection experiments can teach us about the WIMP halo. In other words, we 
want to start the (theoretical) exploration of "WIMP astronomy" . In this first study we use a 
time-averaged recoil spectrum, and assume that no directional information exists. We can thus 
only hope to construct the (time-averaged) one-dimensional velocity distribution j\(v), where v 
is the absolute value of the WIMP velocity in the Earth rest frame. Note that our ansatz is quite 
different from that of the recent paper [22], which assumes a WIMP velocity distribution and 

1 The two experiments might be compatible in some exotic scenarios. One possibility is to postulate rather 
light, m x < 10 GeV, and fast WIMPs with large scattering cross section [20]. Another possible way out is to 
postulate that the detected events are actually inelastic, leading to the production of a second particle that is 
almost, but not exactly, degenerate with the WIMP [21]. 



then analyses with which precision the WIMP mass can be determined from the direct detection 
experiment. 

The remainder of this article is organized as follows. In Sec. 2 we show how to find the velocity 
distribution of WIMPs from the functional form of the recoil spectrum; our assumption here is 
that this functional form has been determined by fitting the data of some (future) experiment (s). 
We then derive formulae for moments of the velocity distribution function, such as the mean 
velocity and the velocity dispersion of WIMPs, which can be compared with model predictions. 
We also discuss some simple halo models. In Sec. 3 we will develop a method that allows to 
reconstruct the WIMP velocity distribution function directly from recorded signal events. This 
allows statistically meaningful tests of predicted distribution functions. We will also show how 
to calculate the moments of the velocity distribution directly from these data. In Sec. 4 we 
conclude our work and discuss some further projects. Some technical details of our calculations 
are given in the Appendices. 



2 One— dimensional velocity distribution function 

In this section we first show how to reconstruct (moments of) the WIMP velocity distribution, 
and then discuss some simple model distributions. 



2.1 Reconstructing the velocity distribution function 

The differential rate for elastic WIMP-nucleus scattering is given by [5]: 
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Here R is the direct detection event rate, i.e., the number of events per unit time and unit mass 
of detector material, Q is the energy deposited in the detector, F(Q) is the elastic nuclear form 
factor, and f\ (v) is the one-dimensional velocity distribution function of the WIMPs impinging 
on the detector. The constant coefficient A is defined as 

A=^, (2) 

where p is the WIMP density near the Earth and a is the total cross section ignoring the form 
fact suppression. The reduced mass m r is defined by 

m r = ^^, (3) 
m x + wn 

where m x is the WIMP mass and m N that of the target nucleus. Finally, v min is the minimal 
incoming velocity of incident WIMPs that can deposit the energy Q in the detector: 



tW = ayJQ, (4) 
where we define 



In Eqs.(l)-(5) we have assumed that the detector essentially only consists of nuclei of a single 
isotope. If the detector contains several different nuclei (e.g. Nal as in the DAMA detector), 



the right-hand side (rhs) of Eq.(l) has to be replaced by a sum of terms, each term describing 
the contribution of one isotope. For simplicity, in the remainder of this article we will focus on 
mono-isotopic detectors. 

In this section we wish to invert Eq.(l), i.e., we want to find an expression for the one- 
dimensional velocity distribution function fi(v) for given (as yet only hypothetical) measured 
recoil spectrum dR/dQ. To that end, we first define 



dF 1 (v) h(v) 



dv v 

i.e., Fi(v) is the primitive of f\{v)/v. Eq.(l) can then be rewritten as 
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Since WIMPs in today's Universe are quite slow, fi(v) must vanish as v approaches infinity, 
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This means that Fi(v — > oo) approaches a finite value. Differentiating both sides of Eq.(7) with 
respect to v min and using Eq.(4), we find 
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Since this expression holds for arbitrary v min , we can write down the following result directly: 
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The rhs of Eq.(10) depends on the as yet unknown constant A. Recall, however, that f\ is 
the normalized velocity distribution, i.e., it is defined to satisfy 

/ f 1 (v)dv = l. (11) 
Jo 

Therefore, the normalized one-dimensional velocity distribution function can be written as 
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with normalization constant Af (see Appendix A) 
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Note that the integral on the rhs of Eq.(13) starts at Q — 0. 



In the next step we wish to compute moments of the distribution function fx. 
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In Eq.(14) we have introduced a threshold energy Qthre- This is needed experimentally, since at 
very low recoil energies, the signal is swamped by electronic noise. Moreover, we will later meet 
expressions that (formally) diverge as Q — > 0. f m in(Qthre) is calculated as in Eq.(4). Inserting 
Eq.(12) into Eq.(14) and integrating by parts, we find (see Appendix A) 
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Physically, (v) is the average WIMP velocity, while (v 2 ) gives the velocity dispersion. 2 We 
emphasize that Eqs.(15) and (16) can be evaluated directly once the recoil spectrum is known; 
knowledge of the functional form of fi(v) is not required. 

Note that the first term in Eq.(15) vanishes for n > if Qthre ~~ * 0. In the same limit, 
(v°) — > J\faI (0)/2 — ► 1 by virtue of Eq.(13). On the other hand, as written in Eqs.(12) and 
(13), the velocity distribution integrated over the experimentally accessible range of WIMP 
velocities gives a value smaller than unity. Using only quantities that can be measured in the 
presence of a nonvanishing energy threshold Qthre , we can replace Eq.(13) by 
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Using A/"(Qthre) in Eq.(12) ensures that the velocity distribution integrated over v > f m in(Qthre) 
gives unity. 

We emphasize again that our final results in Eqs.(12) and (15) are independent of the as 
yet unknown quantity A defined in Eq.(2). They do, however, depend on the WIMP mass m x 
through the coefficient a defined in Eq.(5). This mass can be extracted from a single recoil 
spectrum only if one makes some assumptions about the velocity distribution fi(v). In the kind 
of model-independent analysis pursued here, m x has to be determined by requiring that the recoil 
spectra using two (or more) different target nuclei lead to the same distribution fi(v) through 
Eq.(12). Note that this can be done independent of the detailed particle physics model, which 
determine the value of <To for the two targets. However, one will need to know both form factors, 
which strongly depend on whether spin-dependent or spin-independent interactions dominate 
[5]. Within a given particle physics model, the best determination of m x should eventually come 
from experiments at high-energy colliders. 



2.2 Simple model distributions 

The simplest semi-realistic model halo is a Maxwellian halo. The normalized one-dimensional 
velocity distribution function in the rest frame of our galaxy is then [5] 



2 The dispersion of the function /i in the statistical sense is given by (v 2 ) — (v) 2 . 



where v 220 km/s is the orbital speed of the Sun around the galactic center, which charac- 
terizes the velocity of all virialized objects in the solar vicinity. On the other hand, when we 
take into account the orbital motion of the solar system around the galaxy, as well as the orbit 
of the Earth around the Sun, the distribution function must be modified to [5] 
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where t p ~ June 2nd is the date on which the velocity of the Earth relative to the WIMP halo 
is maximal [8]. Eq.(20) includes the effect of the rotation of the Earth around the Sun (second 
term), but does not allow for the possibility that the halo itself might rotate. 

Substituting these two expressions into Eq.(l), one easily obtains the corresponding recoil 
energy spectra that WIMP direct detection would measure [5]: 
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For future reference we note that the (unrealistic) "reduced" spectrum (i.e., the recoil spectrum 
divided by the squared nuclear form factor) in Eq.(23) is exactly exponential; this remains 
approximately true for the potentially quasi-realistic spectrum in Eq.(24) as well. 

In order to test our formulae, we calculated (v) and (v 2 ), first from the normalized distribution 
functions in Eqs.(18) and (19), 
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Then we used the spectra of Eqs.(23) and (24) in Eqs.(12), (13) and (15), taking Qthre = 0. 
They reproduced the normalized velocity distribution functions in Eqs.(18) and (19), as well as 
the results in Eqs.(25a) to (26b). 

3 Application to experiment 

In the previous section we have derived formulae for the normalized one-dimensional velocity 
distribution function of WIMPs, fi(v), and for its moments (v n ), from the recoil-energy spec- 
trum, dR/dQ. In order to use these expressions, one would need a functional form for dR/dQ. 
In practice this might result from a fit to experimental data. Note that our expression for fi(v) 
in Eq.(12) requires knowledge not only of dR/dQ, but also of its derivative with respect to Q, 
i.e., we need to know both the spectrum and its slope. This will complicate the error analysis 
considerably, if dR/dQ is the result of a fit. 

In this section we therefore go one step further, and derive expressions that allow to recon- 
struct fi(v) and its moments directly from the data. The assumption we have to make is that 
the sample to be analyzed only contains signal events, i.e., is free of background. This should 
be possible in principle: The most copious backgrounds (from radioactive (3 and 7 decays) can 
be discriminated on an event-by-event basis in many modern WIMP detectors. The remaining 
background is dominated by elastic scattering of fast neutrons. While these look just like signal 
events, this background can - at least in principle - be made arbitrarily small by careful shield- 
ing and the use of a muon veto system (to veto cosmic ray induced neutrons). Having a sample 
of pure signal events, we can proceed with a complete statistical analysis of the precision with 
which we can reconstruct fi(v) and its moments. 

In the absence of a true experimental sample of this kind, we had to resort to Monte Carlo 
experiments. To that end we wrote an unweighted event generator. To do so, we had to specify 
the form factor F(Q) appearing in Eq.(l); this is the topic of the first subsection. We then 
proceed to analyze the reconstruction of fi(v) and its moments in two subsequent subsections. 

3.1 The elastic form factor 

We start by presenting the two most commonly used parameterizations of the squared nuclear 
form factor F 2 (Q) appearing in Eq.(l). We focus on spin-independent scattering, which usually 
dominates the event rate; moreover, spin-dependent form factors are still only poorly under- 
stood. 

The simplest form factor is the exponential one, first introduced by Ahlen et al. [23] and 
Freese et al. [8]: 
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where Q is the recoil energy transferred from the incident WIMP to the target nucleus, 



is the nuclear coherence energy and 
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is the radius of the nucleus. 

The exponential form factor implies that the radial density profile of the nucleus has a 
Gaussian form. This Gaussian density profile is simple, but not very realistic. Engel has therefore 
suggested a more accurate form factor [24] , inspired by the Woods-Saxon nuclear density profile, 
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Here j 1 (x) is a spherical Bessel function, 

q = ^2m N Q 
is the transferred 3-momentum, and 
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with 
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In our simulation we used the more accurate Woods-Saxon form factor in Eq.(29). 



3.2 Reconstructing fi(v) 

Since we assume a detector without directional sensitivity, a single event is uniquely character- 
ized by the measured recoil energy Q. Existing experiments such as CRESST [17] and CDMS [16] 
can determine the recoil energy quite accurately. We will see shortly that the statistical errors on 
the reconstructed velocity distribution fi(v) will be quite sizable even for next-generation exper- 
iments, given existing bounds on the scattering rate. We can therefore to good approximation 
ignore the error of Q in our analysis. 

In the following we do not distinguish between the recoil spectrum dR/dQ and the actual 
differential counting rate dN / dQ. Since dRj dQ is usually given per unit detector weight and unit 
time, the two quantities differ by a multiplicative constant. This constant cancels in Eq.(12), 
since it will also appear in the normalization constant (13). 

We divide the total energy range into B bins with central points Q n and widths b n , n = 
1,2, • • • , B. In each bin, N n events will be recorded. Our simulated data set can therefore be 
described by 

Qn- b -f< Qn,i <Qn+ b -f, 2 = 1,2, • • • , N n , 71 = 1, 2, • • • , B. (32) 

The standard estimate for dR/dQ at Q = Q n is then 3 
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3 As usual in the physics literature, we use the same symbol for the estimator, i.e., the "measurement", of 
dR/dQ and for the ideal recoil spectrum. 



The squared statistical error on dR/dQ is accordingly 
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As noted earlier, we also need the slope of the recoil spectrum to reconstruct the velocity 
distribution; see Eq.(12). A rather crude estimator of this slope is 



si, 



d fdR\ 
dQ [dQ 



Q = Qn 



N n {Q > Q n ) - N n {Q < Q n ) 
(fen/2) 2 



(35) 



where N n (Q > Q n ) and N n (Q < Q n ) are the numbers of events in bin n which have measured 
recoil energy Q larger and smaller than Q n , respectively. This estimator is rather crude, since 
it only uses the information in which half of its bin a given event falls. 

It is clear intuitively that an estimator that makes use of the exact Q— value of each event 
should be better. This can e.g. be obtained from the average Q— value in a given bin, defined as 
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Taylor-expanding dR/dQ around Q = Q n , keeping terms up to linear order, gives 
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Eq.(37) allows to calculate Q n as expectation value of Q in the n— th bin: 
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A simple calculation shows that the estimator (39) indeed has a smaller statistical error than 
the one in Eq.(35). The definition (35) immediately implies 
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where we have used Eqs.(33) and (34). In order to calculate the statistical error of the estimator 
S2,n, w e first have to compute 
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Treating the number of events and the average Q— value in a given bin as independent variables 
and using 4 
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This is smaller than the error (40), by a factor 3/4. 

An important observation is that the statistical error of both estimators of the slope of the 
recoil spectrum scale like the bin width to the power —1.5. This can intuitively be understood 
from the argument that the variation of dR/dQ, which we are trying to determine, will be larger 
for larger bins. One would therefore naively conclude that the errors of the estimated slopes 
would be minimized by choosing very large bins. 

However, both Eq.(35) and Eq.(39) reproduce the actual slope at Q — Q n only if the Taylor 
expansion (37) holds; in the presence of higher powers of Q — Q n neither of these estimates 
exactly reproduces the true slope at Q = Q n . Not surprisingly, the influence of these higher 
powers, which induce uncontrolled systematic errors, will increase quickly with increasing bin 
width b n . 

We had seen at the end of Sec. 2 that the predicted recoil spectrum resembles a falling 
exponential. This is confirmed in Fig. 1, which shows the predicted recoil spectrum of a 100 
GeV WIMP on 76 Ge, using the Woods-Saxon form factor (29) and the "shifted Maxwellian" 
velocity distribution of Eq.(19); as usual, we cut the velocity distribution off at a velocity v esc , 
here taken to be 700 km/s, since WIMPs with v > v csc can escape the gravitational well of our 
galaxy. This figure also shows the result of a simulated experiment, where the exposure time 
and cross section are chosen such that the expected number of events is 5,000; these have been 
collected in seven bins in recoil energy. 

While an approximately exponential function can be approximated by a linear ansatz, as in 
Eq.(37), only over a narrow range of Q, i.e., for small bin widths, the logarithm of this function 
can be approximated by a linear ansatz for much wider bins. This corresponds to the ansatz 



while k n is the logarithmic slope of the recoil spectrum in the n— th bin. 

Our next task is to find estimators for f n and k n using (simulated) data. Note that for k n ^ 0, 
f n cannot be estimated from the number of events N n in the n— th bin alone. Instead, one has 
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where we have introduced the dimensionless quantities 
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Figure 1: The curve shows the theoretical predicted recoil energy spectrum (dR/dQ) s ^ for a 
shifted Maxwellian WIMP velocity distribution, and the Woods-Saxon form factor. The data 
points with error bars show simulated experimental data produced from this spectrum (5,000 
total events, m x =100 GeV/c 2 , m N =70.6 GeV/c 2 for 76 Ge, v =220 km/s, v e =23l km/s), and 
galactic escape velocity t> esc = 700 km/s. The vertical error bars show the statistical uncertainties 
of the measurements, while the horizontal error bars indicate the bin widths. 

depends on k n . On the other hand, the second, equivalent expression in Eq.(43) still uses the 
quantities r n = N n /b n as normalization. Evidently, they describe the spectrum dR/dQ at the 
shifted points Q s , n - Equivalence of the two expressions in Eq.(43) implies 

1 , / sinh x„ \ , „ 

Qs,n = Qn + Tr^( " • (48) 

The second expression in Eq. (43) thus has the advantage that the prefactor r n can be computed 
more easily than f n ; on the other hand, while the Q n are simply the midpoints of the n— th bin, 
and can thus be chosen at will, the Q s , n are derived quantities; they depend on the logarithmic 
slopes k n , which we haven't determined yet. 

To do so, we again use the average Q— value in the n— th bin. From Eq.(43) we find: 
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Unfortunately this expression cannot be solved analytically for k n . It is, however, straightforward 
to find k n numerically once Q n is known. Alternatively, we can make use of the second moment 
of the recoil spectrum in the n— th bin, defined as 
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Multiplying both sides of Eq.(49) with b n /x n and adding to Eq.(50), we can calculate the loga- 
rithmic slopes as 
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Note that k n determined either from Eq.(49) or from Eq.(51) is independent of the normalization 
r„ or r n . 

In the following we will estimate the logarithmic slopes from Eq.(49), since it simplifies the 
error analysis somewhat; note that the statistical errors of Q n and (Q — Q n ) 2 n are correlated. 
Using standard error propagation, we have 
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Figure 2: The expected statistical error of our estimator for the logarithmic slope based on 
Eq.(49) as function of the bin width. The bin width is given in units of the inverse logarithmic 
slope, and the error is normalized to its value for b n = —l/k n . In this calculation we have 
assumed that the ansatz (43) is exact within the given bin. 



The error on the average energy transfer can be estimated directly from the data, using 
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where now (Q — Q n ) 2 n is estimated from the data, analogously to the experimental definition of 
Q n in Eq.(36). The second factor in Eq.(52) can be calculated straightforwardly from Eq.(49): 
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where we have defined the auxiliary function 
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For given input values r n , k n and b n , Eqs.(52) and (54) also allow to calculate the expected 
statistical error of the estimated k n , using Eq.(50) to calculate the expected error of Q n — Q n . 
The result is shown in Fig. 2. By normalizing the bin width to the inverse slope, and the expected 
error to its value for a given bin width, the result becomes independent of r n , and can in fact 
be used for all combinations of k n and b n . We observe that for small bins, the expected error 
again scales like b' 1 ' 5 , just as the expected errors (40) and (42) of our two estimators for the 
linear slope. If the bin width is significantly larger than the absolute value of the inverse of the 
logarithmic slope, the error decreases even faster with increasing bin width. 5 

This again argues in favor of using large bins. However, we again have to consider systematic 
errors. After all, it is quite unlikely that the (as yet unknown) recoil spectrum dR/dQ exactly 
satisfies our ansatz (43) over an extended range of Q. Rather, this ansatz should be considered 
as the first terms in a Taylor expansion of the logarithm of dR/dQ. In this case the next-order 
term, which contributes c n (Q — Q n ) 2 in the exponent, will already modify Q n . Since we estimate 
k n from the numerical value of Q n using Eq.(49), which is exact only for c n = 0, any non-zero 
c n will introduce some systematic error in our estimate of k n . 

Fortunately much of this error can be absorbed by a simple trick. According to (43) the 
logarithmic slope is constant over the entire bin, i.e., we could use k n extracted from Eq.(49) as 
estimate of the logarithmic slope at any point Q between Q n — b n /2 and Q n + b n /2. Once c n ^ 
the true logarithmic slope will in fact vary with Q. However, one may hope that the expectation 
value of our estimator still reproduces the true logarithmic slope at some value of Q within the 
n— th bin. 

This is illustrated by Fig. 3, which shows various evaluations of the logarithmic slope within 
one bin as function of the quadratic coefficient c. The true logarithmic slope at the center of 
the bin is, of course, still given by k, independent of the correction c. As argued above, the 
expectation value of our estimator, shown by the dashed (red) line, does depend on c. Note, 
however, that our estimator comes quite close to the true logarithmic slope at the shifted value 
Q s defined in Eq.(48), which is shown by the dot-dashed (blue) line; this is true for both signs of 
c. We therefore conclude that we can minimize the leading systematic error by interpreting our 
estimator of k n as logarithmic slope of the recoil spectrum, not at the center of the bin Q n , but 
at the shifted point Q Stn - Note that Q s>n itself depends on k n ; this, however, does not introduce 



5 Given the exponential form of our ansatz (43) one might assume that the statistical error of the estimated 
values of the k n could be minimized by estimating them from the average values of cxp[n(Q — Qn)], for some 
fixed value of k. For sufficiently small |k| this in fact amounts to using the average value of Q, as described in 
the text. Increasing |k| leads to slightly larger expected statistical errors. 



dR/dQ = exp[k(Q-Q ) + c(Q-Q n ) 2 ], bin width = -2/k 
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Figure 3: Illustration of the systematic error of our estimator for the logarithmic slope based on 
Eq.(49). Here we assume that dR/dQ oc exp[k(Q — Q ) +c(Q — Qo) 2 }, i.e., we amend our ansatz 
(43) by adding a quadratic term to the exponent. The horizontal black line shows the input 
value of k, which now describes the true logarithmic slope only at Q = Q Q . The dashed (red) 
line shows the expectation value of our estimator for the logarithmic slope, while the dot-dashed 
(blue) line shows the true logarithmic slope at the shifted point Q s defined in Eq.(48). The result 
holds for k — — 1 and a bin width b = 2. 

any additional error, if we simply interpret Eq.(48) as an - admittedly somewhat complicated - 
prescription for the determination of the Q— values where we wish to estimate the logarithmic 
slope of the recoil spectrum. 

Using large bins has a second, obvious disadvantage: the number of bins scales inversely 
with their size, i.e., by using large bins we'd be able to estimate f\ only at a small number of 
velocities. This can be alleviated by using overlapping bins, or - equivalently - by combining 
several relatively small bins into overlapping "windows". This means that a given data point 
Q n j may well contribute to several different windows, and hence to the measurement of fx at 
several values of v. This can increase the total amount of information about fx since the only 
information we use about the data points in a given window is encoded in the average recoil 
energy in this window. This averaging effectively destroys information. By letting a given data 
point contribute to several overlapping windows, this loss of information can be reduced. 

A final disadvantage of using large bins or windows is that it would lead to a quite large 
minimal value of v where fx can be reconstructed, simply because the central value Qx, and also 



the shifted point Q S: i, of a large first bin would be quite large. This can be again be alleviated 
by first collecting our data in relatively small bins, and then combining varying numbers of bins 
into overlapping windows. In particular, the first window would be identical with the first bin. 

A final consideration concerns the size of the bins. Choosing fixed bin sizes, and therefore 
also mostly fixed window sizes, would lead to errors on the estimated logarithmic slopes, and 
hence also on the estimates of /i, that increase quickly with increasing Q or v. This is due to 
the essentially exponential form of the recoil spectrum, which would lead to a quickly falling 
number of events in equal-sized bins. We found that we get roughly equal errors in all bins if 
we instead take linearly increasing bins. 

These considerations motivate the following set-up for our mock experimental analysis. We 
start by binning the data, as in Eq.(32), where the bin widths satisfy 

b n = b 1 + (n-l)5; (56) 

here the increment 5 satisfies 
2 



B(B-l) 



TT^Qmax — Qmin — Bb\j , (57) 



B being the total number of bins, and Q max ,min being the (kinematical or instrumental) extrema 
of the recoil energy. We then collect up to nw bins into a window, with smaller windows at 
the borders of the range of Q. In the following we use Latin indices i, j, ■ ■ ■ to label bins, and 
Greek indices /i, z/, • • • to label windows; later on we will use Latin indices a, b, • • • to label 
all events in the sample. For 1 < fi < nw the /i— th window simply consists of the first /i bins; 
for nw < /i < B, the /i— th window consists of bins /i — nw + 1, H — nw + 2, • • • , \i\ and for 
B < a < B + nw — 1, the /i— th window consists of last nw — ^ + B bins. This can also be 
described by introducing the indices and 2 M+ which label the first and last bin contributing 
to the /i— th window, with 

V- = \ „ _ n ,i „ n > V+ = i r „ r > (l<n<B + n w -l). (58) 



ft — n w + 1, /^ > n w 

The center of the i — th bin is called Qi, as before. The total number of windows defined through 
Eq.(58) is evidently W = B + n w — 1- 

The basic observables needed for the reconstruction of f\ are then the number of events iVj 
in the i — th bin, as well as the average Q i defined as in Eq.(36). From these one easily calculates 
the number of events per window, 

v+ 

N„=Y t N i (59) 

as well as the averages 
1 * M+ 

One drawback of the use of overlapping windows in the analysis is that the observables 
defined in Eqs.(59) and (60) are all correlated (for nw ^ 1)- The slope in a given window will 
again be calculated as in Eq.(49), with "bin" quantities replaced by "window" quantities. We 



thus need the covariance matrix for the Q^ — Q^, where is the midpoint of the fi— th window; 
it follows directly from the definition (60): 



V+ r 



cov (q„ - Q„ Q u - Q v ) = —— ]T Ny fe - Q t ) + (Q l - Q„) (q, - Q„) 



,(61) 



where <J 2 (Q i — Qi) is defined as in Eq.(53). In Eq.(61) we have assumed fi < v\ the covariance 
matrix is, of course, symmetric. Moreover, the sum is understood to vanish if the two windows 
fi, v do not overlap, i.e., if i M+ < i v -. 

The ansatz (43) is now assumed to hold over an entire window. We again estimate the 
prefactor as 



N„ 



being the width of the fi— th window. This implies 
1 



cov (r M ,rv) 



(62) 



(63) 



where we have again taken fi < v. Finally, the mixed covariance matrix is given by 



cov (r M , Q u - Q u ) = — — E N i (Qi - QJ) 

W fJ. ly u i=i _ 



(64) 



This sub-matrix is not symmetric under the exchange of fi and v. In the definition of i_ and i + 
we therefore have to distinguish two cases: 



If fi < v : i. 
If fi > v : i_ 



1>V— 1 ^+ ^fJ.+ 



(65) 



As before, the sum in Eq.(64) is understood to vanish if > i + . 

The covariance matrices involving our estimators of the logarithmic slopes k^, derived from 
Eq.(49) with n — > fi everywhere, can be calculated in terms of the covariance matrices in Eqs.(61) 
and (64): 



cov(A; M , k u ) = 



cov (Q» - Q^ Qu - Qv) , 



(66) 



where is as in Eq.(46) with n fi, and the function f(x) has been defined in Eq.(55); and 



cov(r M , k v ) = 



kl 



fix,) 



cov (r M , Q u - Q v ) . 



(67) 



We are now ready to put all pieces together to compute the reconstructed velocity distribution 
and its statistical error. Inserting the ansatz (43) with the substitution n — > fi into Eq.(12), one 
finds the reconstructed velocity distribution 



/l,r(*V) = M 



2Q fi 



F 2 {Qs^ n 

Here, Q s ^ is given by Eq.(48) with n — > fi, and 

= a\fQs,n, 



k. 



(68) 



(69) 



see Eq.(4). Finally, the normalization Af defined in Eq.(13) can be estimated directly from the 
data: 



m-^-t - 



(70) 



where the sum runs over all events in the sample. 

Since neighboring windows overlap, the estimates of fi at adjacent values of t> M are correlated. 
This is described by the covariance matrix 



fiA v »)fiA v ») 



-2M- 



cov(r M , r v ) + (2A0 



Q s ,fj,Q s ,i>T pX v 



Qs,pT p 



/l,rW 



cov^rv) + {fi < — ► v) 



cov(/c M , k u ) 



(71) 



The covariance matrices involving the normalized counting rates r M and logarithmic slopes 
have been given in Eqs.(63), (66) and (67). In principle Eq.(71) should also include contributions 
involving the statistical error of our estimator (70) for M . However, we find this error, and its 
correlations with the errors of the r M and k^, to be negligible compared to the errors included in 
Eq.(71). 

We are finally in a position to present some numerical results. We first validate our results 
by presenting x} distributions, defined via 



w 



J2 C ^ l/i,rW - /i«>] \fiA v ») - AK)] • 



(72) 



fX,V 



Here f\, r is our estimate (68) of the velocity distribution, f\ is the true (input) distribution, and 
C is the inverse of the covariance matrix of Eq.(71). We expect xj to be (roughly) distributed 
according to the standard x 2 distribution when the results of sufficiently many simulated exper- 
iments, with sufficiently many events per experiment, are analyzed. 

Figs. 4 show x} distributions for 5,000 simulated experiments, with on average 500 (top) 
and 5,000 (bottom) events per experiments. Note that the actual number of events in a given 
simulated experiment varies according to the Poisson distribution; otherwise one would introduce 
an artificial correlation between the normalized counting rates in different bins. Moreover, 
the number of bins has been fixed a priori in these analyses. The last bin is typically empty, 
and has therefore been ignored in the analysis. This also reduces the number of windows used 
in the analysis by one, i.e., the upper (lower) frame shows results for W = 6 (W = 8). 

The two histograms in each figure differ by the number of terms that have been included in 
the estimate of the covariance matrix for f 1>r . The solid (blue) histograms have been obtained 
by only including the second term in Eq.(71), while the dashed (red) histograms also include the 
other terms, which are due to the statistical errors on the rescaled event numbers r M . We note 
that including these terms on average leads to a slight overestimate of the true error of /i >r , i.e., 
the average of Xj is somewhat smaller than unity. This is partly due to the fact that we have 
ignored the error on the normalization J\f, which is correlated quite strongly with the errors on 
the r M . 

The lower figure demonstrates that for an average of 5,000 events per experiment the dis- 
tribution of x} values becomes quite similar to the well-known x 2 distribution, shown by the 
smooth curve. At least two effects contribute to the difference. First, we heavily relied on Gaus- 
sian error propagation in our estimate (71) of the covariance matrix of the reconstructed velocity 



5,000 experiments with 500 events each, 5 bins, up to 3 windows per bin 




5,000 experiments with 5,000 events each, 7 bins, up to 3 windows per bin 




Figure 4: The distribution of x} defined in Eq.(72) over 5,000 simulated experiments, with 
on average 500 (top) and 5,000 (bottom) events per experiment. The solid (blue) and dashed 
(red) histograms have been obtained by estimating the covariance matrix for the reconstructed 
velocity distribution excluding (including) the statistical errors on the number of events in the 
windows. The smooth (black) curves show the theoretical \ 2 distributions for the appropriate 
numbers of degrees of freedom. Parameters are as in Fig. 1. 



distribution. This is essentially a Taylor expansion, including only the first non-trivial term. 
It therefore becomes exact only in the limit of small errors, i.e., for large numbers of events in 
a given window. Since the recoil spectrum is falling essentially exponentially, this condition is 
practically always violated at least in the last bin(s), see Fig. 1. We discard windows containing 
less than 3 events, but it is clear that this at best alleviates the problem. In the case at hand, 
this evidently results in an overestimate of the true error. Secondly, our estimator (68) for the 
velocity distribution relies on the estimate of the logarithmic slopes k^, which in turn is based on 
Eq.(49). As illustrated in Fig. 3 this estimate of k^ in general has some systematic error, which 
would tend to increase xj- However, this figure also led us to expect small systematic errors if k^ 
is interpreted as estimator of the logarithmic slope at the shifted points Q s ^. Indeed, as stated 
above, the total expression (71) somewhat overestimates the true error even in the lower frame 
of Figs. 4, which assumes a large number of events but uses a rather small number of bins, which 
thus have to be quite large. Had we instead interpreted k^ as estimator of the slope at Q^, the 
average xj would be about 2.9, indicating that the systematic error would have dominated. 

Figs. 4 also show an excess of simulated experiments with rather large values of xj if the 
covariance matrix for f 1>r is estimated based on the errors on the k^ only. This is true also for 
the upper frame, even though in this case the average value of xj is only about 0.93. To be 
conservative, from now on we therefore take the full Eq.(71) as our estimator of the covariance 
matrix of /i >r , leading to average xj — 0-78 (0.94) for the upper (lower) frame of Figs. 4. 

In Figs. 5 we show results for the reconstructed velocity distribution, for "typical" simulated 
experiments with 500 (top) and 5,000 (bottom) events. In the top frame we choose B = 5 bins, 
the first bin having a width b\ — 8 keV, and combine up to three bins into a window. Since the 
last bin is in fact empty, this leaves us with W = 6 windows, i.e., we can determine f\ for six 
discrete values of the WIMP velocity v; recall that these "measurements" of fi are correlated, 
as indicated by the horizontal bars in the figure. In the lower frame we choose B = 10 bins 
with bi = 10 keV, and combine up to four bins into one window. The bins are thus significantly 
smaller than in the upper frame. As a result, the last two bins are now (almost) empty, leaving 
us with W = 11 windows. 

Figs. 5 indicate that one will need at least a few hundred events for a meaningful direct 
reconstruction of /i. Recall that fx is normalized to unity. The overall magnitude of fx is 
therefore essentially fixed by the range of observed Q— values; only the shape of this distribution 
then remains to be determined. One measure of the information content of the reconstructed 
fi t r is therefore the confidence level with which one can exclude a constant f\. 

In Fig. 6 we show one minus this confidence level, i.e., the probability that a reconstructed 
velocity distribution is compatible with a constant. This has been estimated by defining a 
xj variable as in Eq.(72) for the hypothesis f\ = const., and integrating the theoretical x 2 
distribution over the range x 2 > xj- Here the constant has been chosen as l/(v max — i> mm ), 
where f m m,max have been calculated as in Eq.(4) using the largest and smallest recoil energy, 
respectively, that has been measured in a given experiment. Since this probability differs quite 
widely from one (simulated) experiment to the next, we show both the mean and the median 
probability. We see that we'll need at least 200 events if we want to reject the hypothesis of a 
constant f\ at the 90% C. L. (on average). The confidence level then increases very quickly as 
additional events are added; by the time 1,000 events have been accumulated, we can be quite 
sure that a constant f\ can be excluded with high confidence. 

This confidence level, as well as more general measures of the information that can be ex- 
tracted from a given experiment, depend on the choices of B, b\ and, in particular, nw This is 
illustrated in Table 1, which shows results for different combinations of B, b x and riw for 500 
(first five rows) and 5,000 (last six rows) expected events per experiment. Here the mean and 
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Figure 5: The WIMP velocity distribution reconstructed from a "typical" experiment with 500 
(top) and 5,000 (bottom) events. The smooth curves show the input distributions, which are 
based on Eq.(19). The vertical error bars show the square roots of the diagonal entries of the 
covariance matrix (71); the horizontal bars show the size of the window used in deriving the 
given value of f\ >r . The overlap of these horizontal bars thus shows the range over which the 
values of fi >r are correlated. Parameters as in Fig. 1. 



Average over 1 ,000 experiments 




Figure 6: Estimates of the probability that the reconstructed velocity distribution is compatible 
with being constant, as a function of the average number of events per experiment. These results 
are for optimal combinations of B, b\ and %; they have been obtained by averaging over 1,000 
simulated experiments. The solid (black) and dashed (red) curves show the mean and median 
values of the probability. Parameters as in Fig. 1. 

median probabilities are the same as in Fig. 6. 6 In addition we show the mean of the quantity 
E, defined as 

S = 53C AB/ /i, r (v M )/i ir (v I/ ). (73) 

Formally S determines the significance with which the hypothesis fi = can be rejected. Since 
fx is normalized, this hypothesis is unphysical. Nevertheless £ can be regarded as a measure 
of the information content of a set of reconstructed /i )f .(i> M ); in the absence of correlations, it 
becomes the sum over the inverse squares of the relative errors. Note that, in contrast to xji E 
does not have a 1/W factor in front; after all, by adding more windows we also add more values 

at which f\ is determined, which can increase the information content. 

The first four rows, as well as the last four rows, show the effect of varying nw, the maximal 
number of bins that are collected in a window. We see that there is an optimal choice for this 
quantity. Reducing nw leads to loss of information, as indicated by greatly increased values for 



The results of the first row in the Table have been entered in this figure. 
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5,000 


10 


10 


1 


0.02 


6.7- 10~ 4 


48 
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Table 1: This table illustrates how the binning, and in particular the combination of bins into 
"windows" , affects the information that can be gleaned from the reconstructed WIMP velocity 
distribution. The first four columns show the average number of events in a given experiment, 
the number of bins, the size of the first bin in keV, and the number of bins per window. The 
remaining four columns show the mean and median probability that the reconstructed /i is 
compatible with a constant, the mean of the quantity E defined in Eq.(73), and the average xj 
ofEq.(72). 

the probability that f± jT is compatible with j\ being constant as well as reduced values of E. 
On the other hand, making the windows too large introduces too large systematic uncertainties 
in the estimates of the logarithmic slopes k^, which in turn leads to too large average values of 
xj- This is illustrated by rows four and seven, which have large windows due to our choice of a 
large nw (row 4) or a large bi (row 7). 

The table also shows that the choice of b\ has some impact on the confidence level with which 
the hypothesis of a constant fi can be rejected. We saw in Figs. 5 that our input f\ has a broad 
maximum at v ~ 300 km/s. Rejection of the hypothesis of a constant f\ is therefore optimized 
by maximizing the information about the outer reaches of f\. Getting accurate information 
about fi at large velocities is very difficult; this would need a large number of events at large 
Q, where the counting rate is very small. This leaves the region of small WIMP velocity. By 
choosing a large first bin, one greatly reduces the error on / l r in this first bin, which is also the 
first window; this was illustrated in Fig. 2. In fact, for 500 events and n w = 1 one can formally 
maximize the confidence level with which a constant f\ can be rejected by considering only two 
bins, and making the first bin very large; this is shown in the fifth row. Note that this leads 
to an average xj we U below unity, indicating that in spite of the large bins, systematic errors 
are still insignificant. However, we note that in this case our assumption that the error on M is 
negligible is clearly not justified, since Af receives almost its entire contribution from the large 
first bin. By including the error on r\ but ignoring the (strongly correlated!) error on M we 
clearly over-estimate the total statistical error in this case. Recall also that a large first bin 
leads to a large value for the smallest velocity, v\, where f\ is determined. 

Our "figure of merit" E is less dependent on the details of binning, although, as stated earlier, 
it does strongly benefit from combining several bins into windows. We also note that the optimal 
achievable E is essentially proportional to the number of events in the sample. This is expected, 
since E is something like an inverse squared relative error. 



3.3 Determining moments of fi 



We saw in the previous subsection that a direct reconstruction of the WIMP velocity distribution 
fi will only be possible once several hundred elastic nuclear recoil events have been collected. 
This is a tall order, given that not a single such event has so far been detected (barring the 
possible DAMA observation). The basic reason for the large required event sample is that, f\ 
being a normalized distribution, only information on the shape of j\ is meaningful. In order 
to obtain such shape information via direct reconstruction, we have to separate the events into 
several bins or windows. Moreover, each window should contain sufficiently many events to allow 
an estimate of the slope of the recoil spectrum in this window. 

On the other hand, at the end of Sec. 2 we also gave expressions for the moments of j\. 
With the exception of the moment with n = —1, these are entirely inclusive quantities, i.e., each 
moment is sensitive to the entire data set; no binning is required, nor do we need to determine 
any slope (with one possible minor exception; see below). It thus seems reasonable to expect 
that one can obtain meaningful information about these moments with fewer events. 

An independent motivation for the determination of these moments is that they are sensitive 
also to fi at large values of the WIMP velocity v. We saw above that direct reconstruction of fi 
at large v is very difficult, due to the small number of events expected in this region. Moreover, 
a delta-function-like contribution to f\ at the highest velocity, v = v csc is very difficult to detect 
using direct reconstruction; such a contribution is expected in "late infall" models of galaxy 
formation [11]. 

The experimental implementation of Eq.(15) is quite straightforward. For Qthre — 0, the 
normalization M has already been given in Eq.(70). The case of non-vanishing threshold energy 
Qtbxe can be treated straightforwardly, using Eq.(17). To that end we need to estimate the recoil 
spectrum at the threshold energy. One possibility would be to choose an artificially high value 
of Qthre, and simply count the events in a bin centered on Qthrc- However, in this case the events 
with Q < Qthre would be left out of the determination of the moments. We therefore prefer to 
keep Qthrc as small as experimentally possible, and to estimate the counting rate at threshold 
using the ansatz (43). Since we need the recoil spectrum only at this single point, we only 
have to determine the quantities r± and k\ parameterizing dR/dQ in the first bin; this can be 
done as described in the previous subsection, without the need to distinguish between bins and 
"windows" . Introducing the shorthand notation 

rthre = (S) ' (74) 



the resulting error can be written as 

Qthre — Qs.l ~~ ki 
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The squared errors for r\ and k\ are simply the corresponding diagonal entries of the respective 
covariance matrices given in Eqs.(63) and (66). Finally, the definition (48) of Q s% \ implies 

(If) = *-i + (^) cothx '' (76) 

where x\ = b±ki/2 as before and Q\ is the central Q— value in the first bin. It should be noted 
that the first term in Eq.(15) is negligible for all n > 1 if Qthre — 1 keV; however, even for this 
low threshold energy it contributes significantly to the normalization constant Af, as described 



by Eq.(17). Of course, the first term in Eq.(15) always dominates for n = —1. This is not 
surprising, since the very starting point of our analysis, Eq.(l), already shows that the counting 
rate at Qthre is proportional to the "minus first" moment of the velocity distribution. 
The integral appearing in Eq.(15) can be estimated through the sum 

Q(n-l)/2 

Jn = ?i%p (77) 

see Eq.(70). Since all 4 are determined from the same data, they are correlated, with 

cov(/ n ,/ m ) = £^*|___. (78) 

This can e.g. be seen by writing Eq.(77) as a sum over narrow bins, such that the recoil spectrum 
within each bin can be approximated by a constant. Each term in the sum would then have 
to be multiplied with the number of events in this bin; Eq.(78) then follows from standard 
error propagation. Note that, when re-converted into an integral, the expression for cov(4, 4) 
will diverge logarithmically for Qthrc — ► 0; equivalently, the numerical estimate of this entry can 
become very large if the sample contains events with very small Q— values. Our numerical results 
presented below have therefore been obtained with Qthre — 1 keV; many existing experiments in 
fact require significantly larger energy transfers in their definition of a WIMP signal. 

We also need the correlation between the errors on the estimate of the recoil spectrum at 
Q = Qthre and the integrals I n . It is clear that these quantities are correlated, since the former 
is estimated from all events in the first bin, which of course also contribute to the latter. We 
estimate these correlations again using the ansatz (43), which makes the following prediction for 
the contribution of the first bin to the integrals: 



4,i = n 



fQthrc + bl 
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This immediately implies dl n ^/dr\ = 4,i/?"i, and 
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see Eq.(76). Note that the 4,i and 4+2,1 in Eq.(80) are evaluated as in Eq.(77), with the sum 
extending only over events in the first bin. The correlation we're after is then given by 
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(81) 



These ingredients allow us to compute the covariance matrix for our estimates of the moments 
of the velocity distribution: 



cov 



Af 2 

TV 



m >) 



(v n )(v m )cov(I , Jo) + a n+m (n + \){m + l)cov(/„, I m ) 
- a m (m + l)(u n )cov(/o, I m ) - a n (n + l)(Ocov(/ , Q 
+ D n D m a 2 {r thlc ) - (D m (v n ) + D n (v m ))cov (r thre , J ) 

+ a m (m + l)D„cov(r th rc, I m ) + a n (n + l)L> m cov(r t hrc, 4) 



(82) 



Here we have introduced the modified normalization 



(83) 



which exploits the partial cancellation of the a factors between Eqs.(15) and (17), and the 
quantities 
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(84) 



In our numerical simulations we find that Eq.(77) indeed reproduces the "exact" (input) 
values of the I n if one includes sufficiently many events. In this case it does not matter whether 
one considers a single experiment with a large number of events, or averages over many simulated 
experiments with a relatively small number of events. However, in the second case the average 
values of the reconstructed moments do not exactly converge to the input values. In order to 
understand this, consider the simple case Qthre = 0. The moments are then proportional to 
the ratio I n /Io, see Eqs.(15) and (17). The distortion arises because (I n / Iq) ^ (I n )/(I ), where 
the averaging is over many simulated experiments. In Appendix B we show how this can be 
corrected using Taylor expansion to second order. The leading correction terms for small Qthre 
and not very large first bin are 



S(v n ) = a n Nl{ (n + 1) cov(/ , I n ) - X m I n cov(I , I ) 
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the second line is significant only for n — — 1. Note that this correction becomes very small if 
the statistical errors on the /„ as well as on r t h re become small. 

With this correction, the reconstructed (v n ) indeed closely reproduce the input values after 
averaging over sufficiently many experiments, even if the number of events in a given experiment 
is small. However, the numerical analysis revealed a number of additional problems. These can 
be understood from the observation that the I n in Eq.(77), and even more the entries in their 
covariance matrix (78), receive significant contributions from large Q values, where the counting 
rate itself is already very small. 

This is illustrated in Fig. 7. The a;— axis shows the quantity 



R(Q min ) = \^)dQ } 

jQmin \dQ J 



(86) 



divided by the total counting rate R = R(Q m i n = 0). Here, Q max is the kinematic maximum 
of Q for given input parameters, and Q m i n is varied freely between and Q max - The y— axis 
shows analogously the contributions to some I n (lower curves) and to the corresponding diagonal 
elements of the covariance matrix, i.e., the squared errors (upper curves), that come from the 
region Q > Qmin- In the latter case we have converted the sums in Eq.(78) back into integrals. 
The figure shows that the region of Q— values that contributes 99% of the counting rate only 
contributes about 92% to ii, 73% to J3 and 51% to I 5 ; for the given input parameters, this 
corresponds to the region Q < 103 keV. Even worse, this region only contributes about 35% 
to cov(/ 1 ,/ 1 ) and 5% to cov(/3,/ 3 )! This implies that an experiment collecting only a small 
number of events will typically underestimate (v n ) and, especially, its error; the problem will 
become worse with increasing n. On the other hand, as mentioned above, when averaged over 




Figure 7: This figure has been obtained by introducing an artificial lower bound Qmin in the 
counting rate as well as in the definition of the integrals /„. The x— axis shows the counting 
rate for different Q mm , where small R(Q m m) corresponds to large Q mm . The y— axis shows the 
impact of varying Q m i n on Ii (solid, black), J 3 (dotted, red) and I 5 (dashed, blue); the upper 
curve of a given pattern refers to cov(J n , J n ). The values of the parameters are as in Fig. 1. See 
the text for further details. 



sufficiently many experiments, our estimates for the I n do reproduce the true (input) values. 
This implies that occasionally an experiment will greatly overestimate the I n , the problem again 
getting worse for larger n. 

Our numerical analysis also shows that, after averaging over (very) many experiments, 
Eq.(78) reproduces the mean square deviation between our estimated I n and the true (input) 
value. Nevertheless, we just saw that in most cases this error is being underestimated. In order 
to be conservative, we therefore added "the error on the error" to the diagonal entries of the 
covariance matrix; the off-diagonal entries are then scaled up such that the correlation matrix 
remains unaltered. The squared "error on the error" is defined as 

o- 2 (cov(J n ,J n ))=5:^777T- ( 87 ) 



With this modification, the average x 2 ; again averaged over many experiments, is in the vicinity 
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Figure 8: Estimated moments, and their errors, for a "typical" simulated experiment with 100 
events; recall that the errors are strongly correlated. Parameters are as in Fig. 1. 

of unity at least for the first few moments. 7 

Another problem is that the errors of the /„ are very highly correlated. This can also be 
understood from Fig. 7: A single event at high Q will contribute greatly to all moments with 
sufficiently large n. Numerically we find correlations of more than 98% between (v n ) and (v n+1 ) 
for all n > 2; the correlation between (v ) and (v 3 ) still amounts to more than 87%. This implies 
that the higher moments unfortunately add only little to the available information. Worse, 
attempting to include high moments in a y 2 fit often leads to numerical instabilities; recall that 
a covariance matrix containing 100% correlated entries become singular, i.e., can no longer be 
inverted. In practice only the moments with n ^ 3 therefore seem to be useful. 

We are now ready to present some representative numerical results. Fig. 8 shows the first 
10 moments reconstructed with 100 events, using our standard input parameters (see Fig. 1). 
The estimated values of the moments have been divided by the true values. We see that in this 
"typical" example the high moments are indeed underestimated. We also see that the estimated 
relative errors at first increase with increasing n; this reproduces correctly the trend of the actual 
deviation of the estimated moments from the exact values. However, even after adding "the error 
on the error", we find that the relative errors start to decline again for n > 6. This effect is 
probably entirely spurious; recall that the errors are likely to be even more underestimated than 
the moments themselves. Nevertheless we find it encouraging that already with 100 events a 
couple of moments can be determined with errors of about 15%. 

Figs. 9 show a example of the information that might be gleaned from analyses of recon- 
structed moments of the WIMP velocity distribution. Here we attempt to constrain a possible 

7 Note that (cov(/„. /„)) = ((/„ — I„ !r ) 2 ), where I n ,r are the values of our estimators based on Eq.(77) and /„ 
are the true (input) values, in general does not imply that (cov(/ n , /„)/(/„ — I n ,r) 2 ) = 1< Adding the "error of 
the error" to the covariance matrix brings the average of this ratio closer to unity. 
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Figure 9: x 2 contours, calculated from 3 (top) and 4 (bottom) moments, in the plane spanned 
by the normalization of the "late infall" component and the galactic escape velocity, for typical 
experiments with 25 (top) and 100 (bottom) events. In the upper frame the minimal x 2 value 
is close to the input values Ni = 0, t> esc = 700 km/s; in the lower frame, the location of the 
minimal x 2 is indicated by the star. Parameters are as in Fig. 1. 



•late infall" component in fi [11], denned by the ansatz 
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Here /i ;S h is the standard "shifted Gaussian" distribution (19). As before, we have multiplied it 
with a cut-off at v esc . In addition, we introduce a contribution of WIMPs with fixed velocity, 
which we set equal to t> csc ; these WIMPs are just falling into our galaxy. 8 N t and t> esc are our 
free parameters; N s is then chosen such that the total fi is normalized to unity. We then plot 
contours of 5x 2 , defined as the deviation of \ 2 from its minimal value, where y 2 is defined as 



here (v n ) r are the reconstructed moments in our mock experiment, (v n )(v csc , Ni) are the pre- 
dictions for these quantities based on Eq.(88), and V is the inverse of the covariance matrix 



Figs. 9 show a strong degeneracy in the fit. If the galactic escape velocity f csc is kept fixed 
at its input value of 700 km/s, a quite significant upper bound on the normalization Ni of the 
late infall component could be derived already from our simulated experiment with 25 events. 
However, if v esc is kept free, no significant upper bound can be derived even from the simulated 
experiment with 100 events. Note that the two experiments have been simulated independently, 
i.e., the 25 events used for the analysis in the upper frame are not part of the 100 events used 
in the lower frame. The simulated experiment with 100 events was a bit "unlucky" in that the 
input values lie just outside the A% 2 = 1 contour. As a result, the upper bound on Ni for fixed 
Vesc = 700 km/s actually comes out a little worse in this case than in the experiment with only 
25 events. This is in spite of the fact that the larger data sample allowed us to include one more 
moment in the fit. 

Note that, according to the definition (88), all WIMPs in our galactic neighborhood have 
velocity v < v esc . This implies that a lower bound on t> esc can be derived from the highest observed 
Q— value, f csc > cty / Q + , see Eq.(4). Unfortunately for our standard set of input parameters, this 
method on average only yields lower bounds of about 400 (460) km/s for experiments with 25 
(100) events. Even the experiment with 100 events would then still allow 60% of all WIMPs to 
originate from a late infall component; this is to be compared with the theoretical expectation 
A*"; ^ 10%. Finally, we note that for Ni = 0, it will be essentially impossible to derive a 
meaningful upper bound on v esc from these experiments: because the original "shifted Gaussian" 
velocity distribution is already very small at our input value t> csc = 700 km/s, increasing t> esc has 
essentially no effect on the measured recoil spectrum. 

4 Summary and Conclusions 

In this paper, we have developed methods that allow to extract information on the WIMP 
velocity distribution from the recoil energy spectrum dR/dQ measured in elastic WIMP-nucleus 
scattering experiments. In the long term this information can be used to test or constrain models 
of the dark halo of our galaxy; this information would complement the information on the density 
distribution of WIMPs, which can be derived e.g. from measurements of the galactic rotation 
curve. 

8 Strictly speaking this distribution should also be smeared, since the velocity of the Earth relative to the 
galaxy can add or subtract to the infall velocity. However, as long as v csc is significantly larger than v e , this 
smearing should not matter very much. See Ref. [25] for a discussion of the effect of late WIMP infall on the 
recoil spectrum. 
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To this end, in Sec. 2 we derived expressions that allow to directly reconstruct the normalized 
one-dimensional velocity distribution function of WIMPs, fi(v), given an expression (e.g., a fit 
to data) for the recoil spectrum. We have also derived formulae for the moments of f±. All these 
expressions are independent of the as yet unknown WIMP density near the Earth as well as of 
the WIMP-nucleus cross section; the only information about the nature of the WIMP that is 
needed is its mass. 

Furthermore, in Sec. 3 we have developed methods that allow to apply our expressions directly 
to data, without the need to fit the recoil spectrum to a functional form. We found that a good 
variable that allows direct reconstruction of f\ is the average recoil energy in a given bin (or 
"window"; see Sec. 3.2). This average energy is sensitive to the slope of the recoil spectrum, 
which is the quantity we need to reconstruct f±. We carefully analyzed the statistical errors. 
Unfortunately we found that several hundred events will be needed for this method to be able 
to extract meaningful information on j\. This is partly due to the fact that fi(v) is normalized, 
i.e., only the shape of this distribution contains meaningful information, and partly because this 
shape depends on the slope of the recoil spectrum, which is intrinsically difficult to determine. 

We therefore turned to an analysis of the moments of fc. We found that reliably estimating 
higher moments, and in particular estimating their errors, is difficult. The main reason is that 
these higher moments get large contributions from very rare events with large recoil energies. 
Nevertheless we found that, based only on the first two or three moments, some non-trivial 
information can already be extracted from O(20) events. 

The main emphasis of this exploratory study was on the basic expressions as well as on 
their implementation in actual experiments. The models for fi we tried to constrain in our 
applications (a constant in Sec. 3.2, a "late-infall" component with fixed velocity in the Earth 
rest frame in Sec. 3.3) are not physical; nevertheless they illustrate the difficulties one will have 
in extracting information from these experiments that are of interest for the modeling of the 
galactic WIMP distribution. 

Our analysis is based on several simplifying assumptions. First, we ignored all experimental 
systematic uncertainties, as well as the uncertainty on the determination of the recoil energy 
Q. This is probably quite a good approximation, given that we found that we'll have to live 
with quite large statistical uncertainties in the foreseeable future; recall that not a single WIMP 
event has as yet been unambiguously recorded. 

We also assumed that our detector consists of a single isotope. This is quite realistic for the 
current semiconductor (Si or Ge) detectors. The CRESST detector [17] contains three different 
nuclei. However, by simultaneously measuring heat and light, one might be able to tell on an 
event-by-event basis which kind of nucleus has been struck. In this case, our methods can be 
applied straightforwardly to the three separate sub-spectra. 

Our analyses treat each recorded event as signal, i.e., we ignore backgrounds altogether. At 
least after introducing a lower cut Qthrc on the recoil energy, this may in fact not be unrealistic 
for modern detectors, which contain cosmic ray veto and neutron shielding systems. Background 
subtraction should be relatively straightforward when fitting some function to the data, which 
would allow to use the expressions of Sec. 2. It should also be feasible in the method described 
in Sec. 3.2, if its effect on the average Q— values in the bins can be determined; in particular, 
an approximately flat (Q— independent) background would not change the slope of the recoil 
spectrum. Subtracting the background in the determination of the moments as described in 
Sec. 3.3 might be (even) more difficult. 

As noted earlier, we need to know the WIMP mass m x . This is true for any reconstruction 
method based on data taken with a mono-isotopic detector. In this case one can always "re- 
construct" fi(v), for any (assumed) value of m x . Fortunately in well-motivated WIMP models, 



m x can be determined with high accuracy from future collider data. Even in this case one will 
want to check experimentally that the WIMPs seen in Dark Matter detection experiments are 
in fact the same ones produced at colliders. This can be done by using the methods developed 
here on two different data sets, obtained with different detector materials, and demanding con- 
sistent results for (the moments of) f±. The feasibility of such an analysis is currently under 
investigation. 

In our analysis we ignored the annual modulation of the WIMP flux. Again, given the large 
statistical errors expected in the foreseeable future, this is a reasonable first approximation. 
Nevertheless, eventually one will have to extend the formulae and methods developed here 
to allow for an annual modulation. This is fairly straightforward if the background is again 
negligible. On the other hand, new methods may be needed to extract information on f\ in a 
situation where the total counting rate is dominated by background events; this is most likely 
the case for the DAMA data [15], even if they indeed contain a signal, which remains highly 
controversial. 

In summary, we have begun to explore what direct Dark Matter detection experiments can 
teach us about the velocity distribution of Dark Matter particles in our galactic neighborhood. 
Our analyses show that this will require substantial data samples. We hope this will encour- 
age our experimental colleagues to plan future experiments well beyond the stage of "merely" 
detecting Dark Matter. 
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A Normalization constant and moments of fi 

Since 

v = a^Q, 
we have 
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From Eq.(12) and according to the normalization condition in Eq.(ll), we have, 
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where we have used the conditions 
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Eq.(13) immediately followed from Eq.(A3). 

Using Eqs.(Al), (A2) and integration by parts, we can also find the moments of f\, defined 
with a lower cut-off Qthrc on the energy transfer, as follows: 
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This reproduces Eq.(15) in Sec. 2. 
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B Derivation of the correction terms in Eq.(85) 

Starting point is the observation that we wish to compute the ratio of two integrals, 

Gi _ 1 9i{x) dx _ ^ J2i njgijxj) 
G 2 J g 2 (x) dx Y,j njg 2 {xj ) ' 



(A7) 



In the second step the integrals have been discretized, i.e., replaced by sums over bins % with ni 
events per bin. We now write the rii as sum of average value n-i and fluctuation dnf. 
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Introducing the notation 

G a = J2 n i9a(xi) , a = 1, 2, 

and expanding up to second order in the Srii, we have: 
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We now consider the average over many experiments. Of course, 5rii averages to zero, but 
the product SriiSrij averages to n^-, i.e., it is non-zero for i — j. Hence: 

(^r) - - =2 (j2n i g 1 (x i )g 2 (x i ) ) J + ^ ^^2(2^ . (All) 

The sums appearing in the two correction terms also appear in the definition of the covariance 
matrix between G\ and G 2 - Note that we wish to compute the first term on the right-hand side, 
since in our case the estimators for G\ and G 2 indeed average to the correct values. This then 
leads to the final result 

(=2) cov ( G i> G *) ~ (=1) cov ^ ^) ■ ( Al2 ) 
Applying this result to Eqs.(15) and (17) then immediately leads to Eq.(85). 
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